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CN Abstract 

Motivated by the roll-switching behavior observed in rotating Raylcigh-Benard convection, 

\—A we define a Kiippers-Lortz (K-L) state as a volume-preserving flow with periodic roll switching. 

\,J For an individual roll state, the Lagrangian particle trajectories are periodic. In a system with 

roll-switching, the particles can exhibit three-dimensional, chaotic motion. We study a simple 

phenomenological map that models the Lagrangian dynamics in a K-L state. When the roll axes 

differ by 120° in the plane of rotation, we show that the phase space is dominated by invariant 

tori if the ratio of switching time to roll turnover time is small. When this parameter approaches 

I zero these tori limit onto the classical hexagonal convection patterns, and, as it gets large, the 

^ dynamics becomes fully chaotic and well-mixed. For intermediate values, there are interlinked 

toroidal and poloidal structures separated by chaotic regions. We also compute the exit time 

distributions and show that the unbounded chaotic orbits arc normally diffusive. Although 

the map presumes instantaneous switching between roll states, we show that the qualitative 

features of the flow persist when the model has smooth, overlapping time-dependence for the 

^J roll amplitudes (the Busse-Heikes model). 
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1 Introduction 



S^ In this paper, we study the motion of a passive tracer in two simple kinematic models of 

the Kiippers-Lortz instability, a classical problem in rotating convection. Each model 
is composed of sequentially acting rolls whose axes lie in a plane and are rotated by 
120°. In the limit when the switching from roll to roll is rapid, the tracer trajectories 
limit onto a classical hexagonal convection pattern. As the switching time increases, 
tracer trajectories predominately lie on two- dimensional tori. Eventually, these tori 
are destroyed and the dynamics becomes chaotic and strongly mixing. 

The mixing of tracer elements in a flow is due to a combination of stirring and diffusion. While 
stirring rapidly transports the tracers by kinematic advection, the associated stretching and folding 

*PM was supported by NSF VIGRE grant DMS-9810751. KJ was supported in part by NSF grant OCE-0137347. 
JDM was supported in part by NSF grant DMS-0707659. 



of material lines and surfaces ultimately triggers diffusive processes that homogenize the tracers 
into a blended mixture. Despite this relatively simple description, a complete understanding of 
the spatial and temporal evolution of tracer elements in a given flow remains a mathematically 
challenging problem pQ. This is further compounded if the tracers are not passively advected, 
but nonlinear ly interact with fluid flow and perhaps other species. The phenomena of mixing is 
fundamentally important in many physical and engineering applications and it occurs at a variety 
of scales ranging from the very small (micrometer scale) to the very large (planetary scales and 
beyond). For instance, mixing in microchannels can be used to efficiently homogenize reagents in 
chemical reactions even when the flow is laminar [2]. In reaction-diffusion and combustion systems, 
the efficiency of mixing reactive agents is a crucial factor. Understanding transport for planetary 
scale flows is critical for climate modeling and pollution dispersion in atmospheric science [3] and 
eddy dynamics in oceanography [3] . 

While mixing is prominent in turbulent flows, effective mixing also occurs in laminar flows that 
have chaotic Lagrangian trajectories. The innovative work of Aref [5] gave birth to the subject of 
chaotic advection, which provides a mathematical foundation for investigating mixing and transport 
from a dynamical systems perspective [6] . Aref called his seminal model the "blinking vortex;" it was 
specifically designed to yield nonintegrable Lagrangian trajectories. This model can be interpreted 
as an idealized mixing protocol where passive tracers are successively captured by the velocity fields 
of vortical stirrers that are the analogues of turbulent eddies with finite lifetimes. 

Recently, several models for three-dimensional Lagrangian dynamics have been constructed 
out of vortical roll arrays. Transport in a set of orthogonal, time independent roll arrays was 
investigated in [7J [8]. In this system long-range, superdiffusive transport only occurs along the 
direction perpendicular to the axes of both rolls. This flow could be realized using a combination 
of thermal convection and magnetohydrodynamic forcing. 

In a previous paper, we studied motion in a set of rolls whose axes are orthogonal, and whose 
amplitudes oscillate [9]. This models the oscillations observed in experiments on convection in 
binary fluids [TO] . The experiments showed that the convective instability gives rise to a sequence of 
temporally alternating, orthogonal roll arrays whose axes are parallel to the square boundaries. The 
transition from one set of rolls to the other was observed to be rapid, so that much of the time the 
fluid is dominated by a single set of parallel rolls. In our model, regardless of the choice of temporal 
dependence of the roll amplitudes, the Lagrangian dynamics is restricted two-dimensional surfaces; 
consequently, our model does not lead to three-dimensional mixing, though the two-dimensional 
dynamics can be chaotic. Our model contrasts with that of [7] [8] in that their rolls are 90° out of 
phase. 

Roll-switching was predicted for rotating convection in thin fluids by Kiippers and Lortz [11 | I12|. 
Here a thin fluid layer of large horizontal extent (i.e. small aspect ratio) is rotated about its vertical 
axis. When the rotation rate is small, a stationary roll pattern forms at a critical Rayleigh number. 
Above a critical rotation rate this pattern becomes unstable to a second set of two-dimensional 
rolls whose axes are rotated with respect to the original roll axes by some "switching angle," 0. 
These rolls grow, causing the first set to decay. The secondary rolls then lose stability to a new set 
of rotated rolls and the process repeats. 

This behavior has been observed in a variety of experiments |13| . Typically the weakly su- 
percritical regime shows large patches of homogeneous rolls that loose stability to rolls with an 
axis rotated by a switching angle that varies with the Taylor and Prandtl numbers. For example 



in water (Prandtl number near seven), O = 63 ± 3° [14J. Theoretical investigations confirm the 
dependence of on rotation rate, boundary conditions, and Prandtl number |15| IT6] . Though 
the time dependence of the roll switching is not observed to be precisely periodic, there is a mean 
switching frequency that decreases as the Taylor number increases, but increases with the Rayleigh 
number [17] . 

Motivated by these experiments, we call three-dimensional fluid flow consisting of a set of rolls 
with different axes, whose amplitudes oscillate a Kiippers-Lortz state. In this paper we analyze a 
simple model of the Lagrangian trajectories in a Kiippers-Lortz state. We begin in £J2]by construct- 
ing a time-dependent model with a simple spatial roll structure. When the switching is much faster 
than the roll turnover time, it can be idealized as instantaneous and the flow can be viewed as a 
composition of volume-preserving maps corresponding to the action of each individual roll. This 
gives rise to what we call a blinking roll model in £]3| it is analogous to Aref 's blinking vortex model. 
To check the importance of the assumption of instantaneous switching, we also introduce in £]4] a 
continuous model using the Busse-Heikes system of odes for the evolution of the roll amplitudes 
|18l [19]. When the switching angle is 120° or 60° symmetry, as reviewed in §5] allows us to 
restrict our analysis to a fundamental hexagonal cell. Finally, in £|6] computational results are 
discussed and mixing is analyzed using Lyapunov exponents and anomalous diffusion exponents 

Even when some three-dimensional mixing occurs in these models, it is usually not complete 
because there are families of two-dimensional invariant tori that form barriers to the Lagrangian 
trajectories. Such families of tori are prominent in three-dimensional, incompressible flows; they 
play a roll analogous to that of invariant closed curves in two-dimensional, area-preserving mappings 
|21j . Tori have also been experimentally observed, for example, in a cylindrical container with a 
tilted stirrer |22l [23] . These are visualized with sheet lasers and fluorescent dyes and correlate 
closely with corresponding island chains of the Poincare sections of a model flow. Similar tori 
have also been seen for weakly buoyant tracers in a laminar flow [24 . We investigate some of the 
behavior of the families of tori in our system in §[6j 

2 Kiippers-Lortz Roll Model 

For the purpose of making progress analytically, we use the simplest roll model for rotating convec- 
tion: the linearized solutions in the Boussinesq approximation with stress-free boundary conditions 
|25| . These rolls are tilted by an angle r/ due to the Coriolis force. In scaled coordinates with the 
fluid confined to z £ [tt/2,tt/2], the roll velocity field is 

u^ = A (— tan(ry) sin(y) sin(z), — sin(y) sin(z), — cos(y) cos(z)) . (1) 

Here, A is the amplitude of the rolls and 

tanfo) = j^ , (2) 

where r denotes the Taylor number and k c is the critical wave number at the onset of convection. 
For highly viscous flows and infinite Prandtl number, roll switching starts at a Taylor number 
r c = 54.8 with k c = 3.95, so that r\ ~ 65° [J When the thin fluid layer is not rotating, r = r\ = 0, 

1 In this case, the switching angle for the Kiippers-Lortz state is Q = 59.7° |11] , 



and the velocity field simplifies to: 

uo = A (0, — sin(y) sin(z), — cos(y) cos(z)) (3) 

Although a rigid lower boundary together with a free upper surface is more physical (in which 
case Chandrasekhar functions would be the appropriate building block [25]), this choice would 
complicate both the analytical and numerical treatment. Moreover, in the limit r 3> 1, it has been 
shown |26| [27] that the stress- free solutions converge to the rigid ones outside the passive 0(t~2) 
Ekman boundary layers. 

Trajectories of (fij) are confined to two-dimensional contours about axes parallel to e\ (the x- 
axis). A Kiippers-Lortz state can be modeled in a simple way by using the roll solutions ([I]) as the 
building blocks of the flow. Our model is a superposition of rolls with amplitudes Aj{t) whose axes 
are rotated by angles that are multiples of 9: 

AT-l 

u(x, t)=Y, A j+ i(t)R (j@) u v (it! T (j9)x) . (4) 

j=0 

Here R is the rotation about the z-axis by an angle 9: 

/cos(9) -sin(e) 0\ 
R(&) = sin(9) cos(9) ] . (5) 

V 1/ 

Since each individual roll is incompressible, any composition with arbitrary time-dependence is also 
incompressible. 

We will assume that the rolls act sequentially. That is, when the j roll is active, the amplitudes 
of all other rolls are nearly zero, and that as t increases, Aj(t) decreases to zero, and the next roll 
Aj + i(t) activates and so forth. 

If ¥■ = ~ is rational, then it is appropriate to set N = q, since if the j th roll looses stability to 

the (j + l) st roll, then the whole process starts over after q steps. If, however, J| were irrational, 
then the instability would continue to generate rolls at angles j'9, and the the "number" of rolls, 
iV would be effectively infinite. 

The simplest case of Q corresponds to N = 2 and to orthogonal rolls, as we studied in [9]. 
Here we will study the case N = 3 and a switching angle 9 = 2-7r/3. Letting 

R = #(271-/3) 

we then have the velocity field 

u T (x, t) = Ai(t) Ur? (x) + A 2 (t)Ru ri (R T x) + A 3 (t)R T u ri (Rx) ; (6) 

Note that this model corresponds to a switching angle of 120°, not 60°. However, we shall see in 
£J5] that the roll ([I]) has a reflection symmetry through the origin in the xy plane, thus these two 
are equivalent. 

When the amplitudes are constant and equal, this velocity field gives the (rotating) classical 
hexagonal convection cell 

u h (x) = Ur? (x) + Ru v (R T x.) + R T u v (Rx) , (7) 

The orbits of this autonomous system are either periodic or heteroclinic connecting equilibria on 
the upper and lower boundaries |25j . 



3 The Blinking Roll Map 

The simplest model of sequential rolls corresponds to roll amplitudes that are periodic step func- 
tions. Motivated by the experiments, as well as by simplicity, we assume that both the strength 
and the the activation time of each roll is the same. In this case we can rescale time, setting At — > t 
so that the roll amplitudes become one. The scaled roll activation time, T, thus represents both 
the strength and period of the flow. In this case the roll amplitudes are 

A(t) = { lj ^T(j-l)T<t<jT, 
J \ 0, otherwise 

for < t < NT and periodically repeating thereafter. 

When the amplitude is constant, the flow for the vector field (pi) can be solved exactly in terms 
of Jacobi elliptic functions. To get the flow for the tilted roll, we use the r] = solutions of |9 and 
note note that in (fij), x = i&Ji(r})y, giving 

/ x + tan(r/) [sgn(y) cos -1 S — y] \ 
$ t (x) = sgn(y) cos" 1 S J , (9) 

\ sin" 1 x / 

where 

cos(y)cn(t)dn(t) + sn(t) sin(z) sin 2 (y) 



5(y,z,t) 

x(y,z,t) = 



1 — cos 2 (y)sn 2 (t) 
s'm(z)cn(t)dn(t) — sn(i) cos(y) cos 2 (2;) 



1 — sin 2 (z)sn 2 (t) 
and the elliptic functions have modulus 



k = yl— sin ycos 2 z. (10) 

By setting t = T in ([9]), we get the time-T map 

F(x) = $ T (x) . (11) 



The effective map for the other rolls can be obtained from (11) by rotation. For example in the 
three-roll case the resulting mapping is 

f = (RoFf. (12) 

The map / corresponds to the flow of (pi) over the time 3T. Note that since the vector field is 
incompressible, the map / is volume preserving. Thus, the determinant of the Jacobian, Df, is 
equal to one. 

Different values of can can be easily accommodated by the map by replacing R with R(@). 
When £- is irrational, we simply repeatedly iterate R(@) o F. 



4 Continuous-Time Models 

One simple model of the temporal dependence of the amplitudes of three rolls is the Busse and 
Heikes model [18], 



(13) 



When a and (3 are positive, this system has a heteroclinic cycle between three states. This model 
could also include evolution of the phases of the rolls, but we assume that these are constant. 



A 1 = A X {1- 


l^il 2 - 
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A 3 = 4,(1 - 


\A 3 \ 2 - 


-a\Ai\ 2 - 


-P\M 2 ) 



A thorough review of the dynamic properties of ( 13 ) is given in [19 . When a > 1 > (3, there 
exists an attracting heteroclinic cycle connecting the saddle equilibria, i.e. the pure roll states: 
(A\, A 2 , A3) = (1,0,0), (0,1,0), and (0,0,1). The presence of the attracting equilibria yields a 
system that spends increasing amounts of time around the unstable roll states. This feature is not 
physical, since experimental observations suggest a characteristic period of oscillation associated 
with a Kiippers-Lortz state. Busse and Heikes offered several ideas for overcoming this problem. 



One proposal involves the introduction of small amplitude noise into (13) designed to simulate the 
missing terms in the true amplitude equations. Noisy perturbations of attracting heteroclinic cycles 
were first studied in [28 . Perturbations of this type prevent the system from spending increasing 
amounts of time about any one of the saddle equilibria [19] . The noisy system has a statistically 
"stable" limit cycle with a well-defined mean period. A second suggestion is to add slow spatial 



variation to the amplitude [29 . Adding terms of form d 2 .Ai to each equation in (13) respectively 
also creates a stable oscillation, although its period is difficult to predict. 

We implement the stochastic method by adding the derivative of small white noise terms, £j(i), 



to (13). Specifically, the £i(i) are Wiener processes with mean 2e and variance e 2 . The mean is 
nonzero to help avoid pushing the amplitudes out of the positive octant q] Time series are obtained 
using a second-order stochastic Runge-Kutta solver |30[ [31] . Noise causes the orbit to be pushed 
away from the unstable single roll solutions, preventing it from spending arbitrarily long durations 
near any of those states. 

Detailed pictures of the asymptotically periodic evolution are shown in Fig. [T] The model has 
three parameters a, (3, and e that can be chosen to control mean period, shape, and symmetry of 
the amplitude waveforms. The curves in Fig. [T^ and c show typical realizations of the amplitude 
evolution; the solid curves correspond to much smaller values of e than the dashed curves q In 
Fig. lib and d, probability distribution functions (PDFs) of the periods are shown. For the first 
case the mean period, r sw i tc h is close to 1.5, while in the second it is close to 3. For each case, a step 
size, At, was chosen so that there were roughly 300 steps per period. The PDFs were computed 
using ~ 10 5 periods of a single realization. In both cases, the probability distribution functions for 
the period appear to be nearly Gaussian. 

Let r r oii denote the typical scaled turnover time of a single roll (when A = 1, r ro n = it 2 , see 

2 Though there is a small chance of getting kicked into the unphysical quadrants when the integrations are close 
to the coordinate planes, in practice, choosing mean 2e appears sufficient to prevent this. 

3 fn the numerical work of 96] the waveforms in Fig. Ila and c will be used with tilt angles r\ — 65° and t) = 0° in 
(|6]) respectively. 




Figure 1: Stochastic Busse-Heikes model with parameters given in Table [Tl (a) Evolution of the 
amplitudes: the solid lines correspond to (case a) in Table HI and the dashed lines to (case b). (b) 
Period PDF for an integration time tf = 1.5(10) 5 . The black curve (case a) has r SW itch = 1.500, 
standard deviation a = .0005, skewness £ = —.004, and kurtosis n = —.022. The gray curve (case 
b) has r SW i tc h = 1.500, a = .0008, £ = —.017, and n = .057. (c) Amplitudes for, (case e), solid lines, 
and (case f), dashed lines, (d) Period PDF for tf = 3(10) 5 . For (case e), black, r sw i tc h = 3.00, 
a = .001, C = -057, and k = .002, while for (case f), gray, T switcri = 3.000, a = .003, C = —.006, and 
k = .018. 



Appendix [A]) . The time r SW it c h measures the effective strength of the rolls; consequently when 



''"switch 
Troll 



« 1 



(14) 



the strength of the rolls is relatively weak. As this parameter grows the rolls become stronger. 

Let rtrans denote the typical transition time between roll states (e.g., the transition time between 
0° and 120° states). If the ratio of transition time to switching time is small: 



Ttr 



« 1 



^switch 



(15) 



then a relatively large time is spent in a particular roll state compared to the transition to the next 
roll. The parameter e is the key to controlling this ratio; indeed, comparing the solid and dashed 
curves in in Fig. [T] it is clear that as e grows this inequality becomes invalid. 



5 Geometry 

Since the vertical component of the velocity (fTl) vanishes at z = ±7r/2, the orbits are confined to 
the infinite domain 

(16) 



■D={(x,y,z):--<z<-j . 



As we will review below, the velocity field mty is invariant under a discrete translational sym- 
metry generated by the unit vectors 



«i = (1,0,0) 



U 2 



1 v^ 

2'~2~ 
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y= 271 
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= -271 



Figure 2: Schematic of the roll geometry in z = plane. The open circles correspond to saddle 



equilibria at z 



The solid lines correspond to roll boundaries and the dotted lines correspond 



to roll axes. The + signs denote right-handed or positive rotation about its axis. The hexagon, H 
is a fundamental cell of the flow. 



These vectors generate the planar "wallpaper" group denoted by p3: 



p3 = < t mn (x) = x + —={mui + nu 2 ) : m,n eZ 



(17) 



corresponding to a single, order-three axis of symmetry. Here we have scaled the length of the 
vectors to correspond with the unit cell size in our examples. This group has a fundamental cell 
that is a regular hexagon, TC, of width 4-7r/\/3 and height 87r/3, see Fig. lA 

Regardless of the choice of time dependence for the amplitudes, the vector field (pi has a fourteen 
saddle equilibria on the upper and lower boundaries, z = ±|; these are shown as the dots in Fig. 12} 
Each vertically separated pair is connected by an invariant line giving a heteroclinic connection. 
For example, the saddles at (0,0,±|) are connected by the invariant line {x = y = 0}. Along 
this heteroclinic connection, the fluid motion is downward. The remaining twelve saddles can be 
found on the boundary of TC, six at z = tt/2 and six at z = —tt/2. There are heteroclinic orbits 
connecting these pairs that move upward. For example, there is an upward moving heteroclinic 



connection between the points ( -j= , 



to 



2tt 
v/3 



,o, 



Note that since the map ( 12 ) is an exact solution of <\6h for the special case of step function am- 



plitudes, the equilibria and heteroclinic connections of the flow become fixed points and heteroclinic 
lines of the map. 



6 Numerical Explorations 

In this section we study the mixing and transport for the flow (I6| and compare it to that for the 



blinking roll mapping (12). One local measure of mixing corresponds to the Lyapunov exponents 
as a function of initial condition. There is a correlation between large exponents and regions of 
strong mixing. On the other hand in regions with many invariant tori, the exponents are zero or 
near zero. 

A second measure is the rate of particle drift throughout the domain. A first goal is to distinguish 
between bounded and unbounded orbits. We will study those bounded orbits that remain within 
the fundamental domain TC — we call these orbits "trapped." Generally the "regular" orbits appear 
to be trapped (for example, those on invariant tori), especially for the case T\ = T2 = T3; however, 
there appear to be bounded invariant structures larger than a single fundamental domain when the 
mapping times are not equal. We will also compute the exit time distribution from a fundamental 
domain. Transport of unbounded orbits can be quantified by computing the growth rate of distance 
with t and comparing this to diffusive behavior. Since the trajectories are bounded in z, we measure 
the drift in (x, y) by looking at the growth rate of r 2 = x 2 + y 2 . 

6.1 Lyapunov Phase Portraits 

In this section, we compute the three Lyapunov exponents, (/xi, (12,^3), as a function of initial 
condition for both the flow and map models. Ideally we would compute the exponents for each 
initial condition on a three-dimensional mesh covering the fundamental cell, but this would be 
prohibitively expensive in computation time and the results would be difficult to display. Instead, 
we use two, two-dimensional slices in the fundamental cell TC: the horizontal midplane, z = 0, and 
the vertical cross section, y = 0. The horizontal slice is hexagonal and is divided into a 70 x 70 
grid and the vertical slice is rectangular and divided into a 100 x 50 grid. 

Since the map is volume-preserving and the flow is incompressible, the sum \i\ + \i% + ^ = for 
both models. In Fig. |3j we show only the value of the largest exponent /ii for each initial condition 
on these grids. We call these pictures "Lyapunov phase portraits;" they help to differentiate between 
chaotic and regular regions and provide a local measure of chaos that correlates with rapid mixing. 
The exponents are computed using the QR method of Rangarajan et al |34[ 135] . The second and 
third Lyapunov exponents were also computed for all cases. Typically we observed that I//2I ~ 10~ 4 
which is two orders of magnitude smaller than fj,i in the chaotic region. Given that only 10 5 
iterations were used in the computation and that converge is slow, a longer integration might be 
expected to yield an evan smaller second exponent. Thus our results are consistent with \i<i = 0, 
and consequently with ^3 = —[i\. Note that an autonomous flow always has a zero exponent. More 
generally, if the flow or map were reversible then we would also expect /i2 = 0. We have not been 
able to show, however, that our systems are reversible, and thus do not have an explanation for 
the near vanishing of \i% . 

By computing the Lyapunov phase portraits for both the flow ([6]) and the map (12), a direct 



comparison can be made between the two dynamical systems if the parameter values are equivalent. 
For the flow we use a full integration of the stochastic Busse-Heikes model to generate the amplitude 
waveforms; those shown in Fig. [I] give a typical cycle. Since there are three blinking-rolls per map 
iteration, the corresponding parameter T for the mapping is one third of the average period of the 
amplitude functions. 



Table 1 : Parameters for the Lyapunov phase portraits of Fig. [3] 
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The results are shown in Fig. [3] and the corresponding parameter values are given in Table [T] In 
each case the total time corresponds to approximately 10 5 oscillation cycles. The average period of 
oscillation is approximately 1.5 for the first two flow cases (panels (a) and (b)), and we choose the 
corresponding mapping time in panel (c) to be 0.5 so that the total integration time represented 
by iterating the map 10 5 times is the same as that for the flow. For panels (e) and (f) the average 
period is 3.0, and the corresponding mapping time T = 1.0 is used in panel (g). f\ 

Results for the flow with r\ = 65° are displayed in Fig. [3^, and Fig. [3J3. The amplitude waveforms 
in these two cases differ in their shapes, as shown by the solid and dashed curves in Fig. [T] The 
images indicate that there are toroidal regions of near zero Lyapunov exponent surrounded by a 
chaotic sea. The breaks in the toroidal region on the hexagonal slice correspond to locations where 
the torus does not intersect this slice. It is striking that the tori exist even when the three amplitude 
waveforms have significant intervals of overlap as they do for Fig. [3]x This suggests that the details 
of the amplitude wave form is not a critical factor in the existence of the tori. The corresponding 
Lyapunov phase portrait for the map, shown in Fig. pp, is visually quite similar to (and as shown 
in Table^lhas a high correlation with) the first two flow cases. In Fig.^e and Fig.^r, the Lyapunov 
phase portraits are shown for the waveforms in Fig. [TJ: with r] = 0°. Once again, these portraits 
show a high correlation with those for the corresponding mapping, shown in Fig. |3^. Although 
the T) = 0° case is nonphysical for rotating convection, it is interesting to see that the toroidal 
structures persist for small rj and are robust as this parameter varies. 

A few orbits inside the regions of low chaos (red, blue, and green) and one in the chaotic region 
(purple) are shown in Fig. [3h for (12). The regions of near-zero Lyapunov exponent correspond 
to a set of nested two-dimensional tori. At the center of this nested set is apparently an elliptic, 
invariant curve. The green orbit in Fig. [3]i appears to correspond to a torus that encloses an 
invariant curve that is helically wound about the central family of tori. 

In each case, the close correlation between the Lyapunov exponents of the flow and mapping 
is clear from the pictures. To quantify this connection, we computed the correlation coefficient 
between the values of ji\ for mapping and for the flow on each slice. The results are given in 
Table [2J coefficients near one indicate similar structure whereas coefficients near zero suggest large 
differences between the data sets. In each case, a reasonably strong correlation is obtained. It 



4 The average computation time was about 4 cpu days for the flow. Computations were performed on a 4 processor, 
64 bit, Itanium 2 1.3GHz (32KB LI cache, 256 KB L2 cache, 3MB L3 cache) architecture. In contrast, the Lyapunov 
exponents for the discrete maps took 2 — 3 cpu hours on the same machine. 
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Figure 3: Lyapunov phase portraits for the flow (pi) and map (12) for the parameter values given 
in Table [X] The largest Lyapunov exponent is computed by integrating or iterating for 10 5 periods 
and the graphed values are obtained by averaging over the final 100 periods. The color scale for 
[i\ is indicated at the bottom of each pane. Sample trajectories of the mapping are given in panels 
(d) and (h) showing regions occupied by invariant tori. Panels (a),(b), and (c) have the same color 
scale as do panels (e), (f), and (g). 



Table 2: Correlation Coefficient of largest Lyapunov exponent between the map and the flow 
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is notable that as the ratio (15) increases, the correlation drops PI This seems to predominantly 
occur near the outer boundary of the regular regions where the step function approximation to 
the amplitudes does a poor job of approximating the flow: the map tends to have a larger region 
foliated by invariant tori than does the flow. 

As computations with the map are much quicker than those for the flow it is easier to investigate 
the effect of changing parameters on the map. Moreover, changing the switching angle from 



120° is problematic for the flow since the three-mode model (13) no longer applies. It is very 



easy, however, to change the switching angle in the map by choosing the rotation angle for the 
transformation R in ( |12[ ) . One of the immediate observations is that tori are robust features of the 
mapping as the switching and tilt angles and time step T are varied. For example, a large family of 
tori can be seen Fig. Hk (viewed from a camera looking down from large z). For this case the ratio 
( 14 ) is small. The size of the innermost torus shown is on the order of one hexagonal convection 



cell; however, since O ^ 120° the system no longer has hexagonal symmetry. Outside these tori, 
the motion is chaotic and diffusive (not shown in the figure). 

A final example in Fig. [4J3 shows a number of helical tori enclosing the central family. In the 
center of the light gray orbits are invariant curves that wind twice in the toroidal direction before 
returning to their original positions. It is possible that these families of tori are created by period- 
doubling bifurcations of the invariant curve at the center of the main torus or a resonant bifurcation 
of a two-torus, but we have not verified this. 

There are too many parameters in the system to completely explore the full range of possible 
dynamics. A critical point is that the discrete symmetry does not produce the tori: these structures 
are features of the flow for many parameter values. 

6.2 Breakdown of Integrability 

When the time T is very small, it is easy to see that the map ( |12[ ) is, to first order in T, the same 
as the time-T map of the autonomous field ttfh. Thus for a fixed total iteration time, we expect 
that the map dynamics will limit on the integrable flow of (17b as T — > 0. Indeed, we observe that 
the mapping shows the same closed, integrable trajectories of the classical hexagonal convection 
patterns [25] when T approaches zero for a given fixed total time nT. 

As T grows, the closed trajectories are destroyed and are replaced by nested invariant tori 
separated by small zones of chaos as illustrated in Fig. [5] This example shows six families of 

5 Here Ttrans is computed by finding the total time t that Ai(t) > .95 or Ai(t) < .05 for i — 1, 2, 3. This value is 
subtracted from T sw i tc h- 
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Figure 4: Phase portraits of (12) for (a) T = .2, rj = 65°, and 9 = 122°, and (b) T = 1, rj = 0° 



G = 45°. We plot points of the map RF for each iterate instead of every third iterate as in (12 ). 

c 






Figure 5: Phase portraits of (12) for = 120°, rj = 0, and T = 0.2. Panel (a) shows the Lyapunov 



phase portrait, (b) the elliptic invariant curves at the center of the poloidal tori and (c) a chaotic 
trajectory in the region between the tori. 



"poloidal" two-tori linked with a central "toroidal" family. The poloidal tori can be seen crossing 
the z = section in the Lyapunov phase portrait of panel (a) near the center and edges of the 
hexagon. Invariant curves at the center of the poloidal and toroidal families of tori, shown in panel 
(b), can be found using the simple algorithm discussed in Appendix |B| A trajectory in the chaotic 
region between the tori is shown in a top down view in panel (c). This domain can also be seen as 
the light grey regions in the Lyapunov phase portrait of panel (a) . 

This same structure is seen when T is small for any value of r/. As rj increases the tori becomes 
increasingly twisted in the same manner as the classical hexagons. For each rj, as T is further 
increased, the six poloidal tori disappear leaving only the toroidal family. Eventually the toroidal 
family is destroyed as well and the cell appears to be completely chaotic (e.g. Fig. [6]i below). In 
Table [3j the largest value of T for which the map has visible tori is given as a function of rj. For these 
T values, the central invariant curves are warped and twisted (see Fig. pk) and the surrounding 
two-tori have a small poloidal diameter. For large rj, the transition from integrability to complete 
chaos is rapid. 
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Table 3: Largest time, T, for which (12) has invariant tori as rj varies 



T] 





5 


15 


25 


35 


45 


55 


65 


75 


85 


T 


1.81 


1.8 


1.7 


1.6 


1.5 


1.35 


.97 


.6 


.6 


.25 
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Figure 6: Toroidal invariant curves of (12) for r\ = 0. (a) Curves for T = 0.2,0.6, 1, 1.4, and 1.8. 



The solid black line is (18). (b) Rotation Number to as a function of mapping time T. (c) Phase 
portrait for T = 2 where there are no invariant tori giving complete mixing. 



The evolution of the toroidal invariant curve at the center of the toroidal family of two-tori is 
shown in more detail in Fig. |6j These curves are found using the center of mass algorithm algorithm 
of Appendix pi In panel (a), the curve is plotted for a number of different parameter values. As 
T — ► the invariant curve approaches the closed curve of the autonomous flow of (|7| given by 



cos(y) + 2 cos 



V^: 



X 



COS 



. 



(18) 



This is shown as the darkest black curve at z = 0. For larger T values the invariant curve becomes 
increasingly distorted, and appears to lose differentiability. The curve is destroyed for T ~ 1.81. 
For larger values of T, as shown in panel (c), the phase portrait appears to be completely chaotic. 
Panel (b) shows the longitudinal rotation number of the invariant curves, u, as a function of T. As 
T — » 0, to — > since the mapping limits to a time T map of the autonomous flow. As T increases 
the rotation number grows monotonically. There are gaps in the rotation number corresponding to 
parameter intervals where the curve is destroyed and reforms. 

The six poloidal tori shown in Fig. [5] also appear to arise from streamlines of the autonomous, 
hexagonal flow ([7]) as T increases from zero. Symmetry implies that they appear in the planes 
x = and y = and their rotations by tt/3. When r\ = 0, the streamlines of (fn) in these planes 
correspond to contours of the functions 



*i(0,y,*) 



cos(z) sin 2 



cos(z) sin 3 




(19) 
(20) 
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respectively. When rj ^ 0, a closed form expression for the streamlines is difficult to find. Using 
T = 10~ 4 and the center of mass algorithm, computations show that the invariant curves at 
the center of the six elliptic tori limit onto the streamline defined by M/i ~ 0.3868 This streamline 
presumably is selected by a mechanism similar to the classical Poincare-Birkhoff island chain created 
in when an integrable area-preserving map is perturbed. 

Just as in the two-dimensional, island chain case, one of the ^2 streamlines apparently becomes 
a hyperbolic (saddle) invariant curve in the same manner that the ^1 = 0.3868 streamline becomes 
an elliptic invariant curve. Indeed, numerical computations show localized chaos around the planes 
y = and y = ±y3x (recall Fig. pa and c) that grow as T increases from zero. Transverse 
intersections of the stable and unstable manifolds of these saddle curves provide a mechanism for 
the onset of three-dimensional chaos [46 . 

The disappearance of the elliptic invariant curves as T increases is expected to be caused by 
various resonant bifurcations between their longitudinal and transverse rotation numbers. This 
hypothesis is supported by the observation that the regular region surrounding the elliptic curves 
shrinks as they approach destruction. We hope to report further on these bifurcations in a future 
paper. 

6.3 Transport 

While Lyapunov exponents measure the divergence of neighboring trajectories, they give no infor- 
mation about how trajectories disperse over macroscopic distances, that is about "transport." The 
goal of this section is to develop several measures of transport. 

A first distinction to make is between trajectories that are bounded and those that are not. 
One class of bounded orbits are those that remain forever trapped in a single fundamental cell, 
TC. These, for example, include some of the invariant tori that we previously studied, though it 
is possible for invariant tori to exist on scales larger than a single cell. Since some trajectories 
never leave TC, it is reasonable to restrict the calculation to the set of initial conditions whose 
trajectories eventually do leave; this set is equivalent, up to a set of zero volume, to the subset that 
is "accessible" to trajectories that begin outside TC [36]. The accessible set can be constructed by 
considering the incoming set TCi n which is defined to be 

H in = H\ F{TC) , 

i.e., the portion of TC that does not intersect F(TC). According to this definition, F (Tii n ) is outside 
of TC, thus every point in TCi n has just entered TC. The accessible portion of TC is the portion that 
can be reached by starting in the incoming set: 



TCacc=\jF\TC m )r\TC 



t=0 

In order to leave TC, a trajectory must first land in the exit set 

TC ex it = H\F-\H). 
By volume preservation, the incoming and exit sets have the same volume. 
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Figure 7: Exit time and mean square displacement with parameters given in Table II] (a) Exit time 
distribution for a 1000 x 1000 x 500 grid in fundamental cell, (b) Mean-squared displacement for 104 
randomly sampled points in TC ex it- (c) Mean-squared displacement for the flow mj with amplitudes 
given in for 103 randomly sampled points in TC ex it- 



A simple measure of transport is the "exit time distribution," ip(t), defined to be the probability 
that a trajectory that starts in TCi n at t = 0, first leaves TC at time t. To compute this we cover TC 
with a grid of points. Each point is iterated with the inverse map to see if it remains in TC — if so, 
it is discarded. The remaining points are in the incoming set, and they are iterated forward until 
they exit TC. 

Exit time distributions for hyperbolic systems are typically exponential; however, it is common 
in volume-preserving dynamics for ip to have a power-law form |21| : 



ip(t) ~ f 



(21) 



For the case of a volume-preserving map, it is a consequence of Kac's theorem that the average exit 
time 

roo 

Texit= / tll>(t)dt, (22) 

JO 

exists |36j . Thus ip(t) = o{t~ 2 ) as t — ► oo. The computations of exit time distributions, shown in 
Fig. [7^, are consistent with (21) and 7 pa 2.1 — 2.5. In each case, as the theory requires, T ex u < 00 
since the distribution has power law tail with 7 > 2. 

Chaotic trajectories that leave the fundamental cell tend to have unbounded orbits (though as we 
have seen there can also be superstructures extending over several copies of TC). One measure of the 
transport properties of unbounded trajectories is the growth rate of the mean-squared displacement, 



i.e. 



(r 2 (t))~t» 



(23) 



The case fj, = 1 corresponds to diffusive motion; however, many dynamical systems exhibit either 
subdiffusive (/j, < 1) or superdiffusive (/x > 1) behavior [37]. Anomalous diffusion has been observed 
in many systems, for example, in two [39114*0] and three-dimensional fluids [7\, plasma physics [38 , 
and electron transport on semiconductor super lattices |41j . 
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One phenomenological model that gives anomalous diffusion of the form of (23) is a continuous 
time random walk. In this model particles are assumed to move from one point to another under 
a sequence of "Levy flights" whose length and duration are taken from probability distributions 
with power law tails. Such models display both normal and superdiffusive behavior |42| I33~| 133] . A 
trapping probability distribution, separating jumps by trapping events in which the particle remains 
at its current position for a randomly selected time can also be added to the model — in our case 
this distribution is ip(t). For this situation, the dynamics can also exhibit subdiffusive behavior. 
Subdiffusion occurs when the motion is dominated by trapping events, and superdiffusion when 
flights are the more influential factor |45j . 

For our model the average exit time is finite; consequently, if there are no flights, one can 
expect normally diffusive behavior |47j . Though it is not easy to directly measure the distribution 
of particles undergoing so-called "Levy flights," we can compute the mean squared growth rate 



fi, (23). To do this we randomly sample 10 points in H e xit, iterating each 10 times. We also 
compute the mean square displacement for the continuous system, though in this case, since much 
computation is involve we only use 10 3 initial points and integrate for 10 4 periods. In both cases 
the results, shown in Fig. [7] are consistent with [x = 1 however \x = 1 does not always fall within 
the 90% confidence interval. However, the least squares technique for computing the asymptotic 
growth rate may not give an accurate estimate of the error. In general though, it appears that these 
models are normally diffusive and that there are no Levy flights, unlike many similar dynamical 
systems. 

7 Discussion 

We have studied two idealized, kinematic models of the Kiippers-Lortz state. The building blocks 
of the models are vortical roll arrays whose axes are rotated through an angle O with respect to one 
another. In the continuous-time construction, the time dependence of the amplitude of the rolls 



is controlled by a stochastic Busse-Heikes system (13). When the amplitudes of the rolls are step 



functions and the rolls switch instantaneously from one set to another, we derived an analytical 



map, (12). 

At the simplest level, the transition to chaos is governed by the the ratio of switching time, 
^switch; to turnover turnover time, T ro ii- As this parameter goes to zero, the flow approaches the 
autonomous, integrable case and when it is sufficiently small, the dynamics is similar, showing the 
classical hexagonal patterns observed in rotating Rayleigh-Benard convection. 

As this ratio grows, both the flow and the map have a family of invariant tori surrounding an 
elliptic invariant curve. Outside these tori is a chaotic sea. The tori undergo a complex sequence of 
bifurcations as the parameters change that we hope to report on in a future paper. In some cases 
we have also observed localized mixing for trajectories trapped within an invariant torus. This 
occurs, for example, if the central invariant curve is temporarily destroyed by a bifurcation, but 
the surrounding two-tori persist. 

For large values of the the switching time, the Lyapunov phase portraits indicate nearly com- 
plete mixing. We characterized the transport properties of our system by computing the mean 
exit-time from the fundamental cell and the mean-squared drift of trajectories over the extended 
domain. Computations for the blinking-roll mapping and the flow showed algebraic decay of the 
exit time distribution, and the mean-squared growth consistent with normal diffusion. These re- 
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suits, in conjunction with anomalous diffusion theory, suggest that there is no analogue of the open 
streamlines in an autonomous flow or the "accelerator modes" of area-preserving maps that would 
effect long range, superdiffusive transport. 

We also showed that the dynamics of the flow and map models is highly correlated even when 
the switching time is far from instantaneous. This observation suggests that a much broader class 
of continuous-time models, perhaps with periodic switching between more than three states, can be 
accurately modeled with the blinking roll mapping. With this in mind, one could easily extend the 
blinking roll system to include more realistic spatial descriptions for the flow, like Chandrasekhar 
functions, which are solutions to the linearized Boussinesq approximation with a no-slip lower 
boundary condition. However, it would be more difficult to obtain an explicit mapping and the 
computations would consequently be slower. For the simple case explored here, the Lyapunov phase 
portraits can be computed about twenty times faster for the map than for the flow. 

The mapping has tori over a wide range of switching and tilt angles. We explored some of their 
behavior, but many mysteries remain. For instance, the middle Lyapunov exponent appears to 
be zero to numerical accuracy for both the continuous and discrete models. While this would be 
rigorously true for autonomous flows or reversible maps, we do not have an explanation of this for 
our systems. In addition, it would be interesting to study in detail the bifurcations that give rise to 
the families of invariant tori and curves, as well as those that ultimately lead to their destruction. 
Numerical studies show that destruction of invariant curves and tori typically enhances the mixing, 
but as of yet there is no satisfactory theory of transport for these three-dimensional systems. 

A Appendix: Mean Roll Turnover Time 

Here, we compute the average roll turnover time, T ro ii, for a single roll with streamfunction ip = 
sin(y) cos(-z). Suppose an initial condition is placed along the line z = 0, < y = yo < §• The flow 
along this streamline is given by (T9J) (with r\ = 0): 



( . . i / cos(yo)cn(t)dn(t) \ _, 

®t {x, 2/o, 0) = [x, sgn(yo) cos ^— — — ^—r , sin (-sn(i) cos(y J) 

V \l-cos z {yo)sn z (t)J 



(24) 



Now, solve for the t value such that 



(7T 7T 

This corresponds to one quarter of a rotation around a given streamfunction. This is possible since 
the tracers rotate counterclockwise in the domain {(x,y) : < y < it, — \ < z < |}. From the third 
component, we find a simple equation for t: 



Here, the modulus of the elliptic functions is (10) and F(%,k) is the complete elliptic integral of 



the first kind. To get the average roll turnover time for all streamfunctions, we integrate over k to 
find 



4jf I F(=,*)«-r'. 



Troll 

It is straightforward to show that T ro ii is independent of rj. 
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B Appendix: Finding Invariant Curves 

Locating invariant curves is a more challenging task than finding fixed or periodic points of a 
map. For a map, the dynamics on an invariant curve is typically conjugate to a rigid rotation with 
irrational winding number. Therefore, one cannot use standard root finding techniques since the 
orbit never returns to itself. Though there are many standard algorithms for finding periodic orbits, 
methods for finding quasiperiodic ones are not as well developed. One algorithm uses Newton's 
method coupled with an interpolation algorithm to find a single point on an invariant curve |49J . 
In another, the curve is approximated by a finite Forier series and Newton's method is used to 
solve for the coefficients [50j . This this method is hindered by the assumption of conjugacy to rigid 
rotation number with an a priori known rotation number. 

We use a simple algorithm that computes a single point on the invariant curve; however, only 
one point is needed since the rest of the invariant curve can be found through iteration. The method 
is limited because it only works for the case of an elliptic invariant curve, and it is rather slow. On 
the other hand, the method has significant advantages in that it does not require foreknowledge of 
the rotation number u and it is easy to implement. 

The first step is to choose an appropriate plane that is expected to intersect the invariant curve 
transversely. We fatten this plane into a thin slice of width e and choose a point xo in the slice close 
to where the curve is expected. Iterate the initial point and record the times it lands within the 
slice. Once N such near returns are found, we compute the center of mass of these points, x. If the 
elliptic curve is surrounded by a family of two-tori, and these have a convex intersection with the 
slice then x will approximate the position of the enclosed invariant curve. We now decrease e and 
repeat the process with xo = x as the new starting position. This procedure appears to work well 
to find a point on the invariant curve to machine precision. Given this point, we can now iterate 
to find the longitudinal rotation number to of the curve. 
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